DIPGEOPR_03

class: inverse,left, middle
background-image: url(background.png)
background-size: cover

##Módulo 3: Percepción remota, análisis masivo y Google Earth Engine.

### Pre-procesamiento satelital I

Gabriel Castro

mail gabriel.castro.b@pucv.cl</a>

.large[<a href="https://www.pucv.cl/uuaa/site/edic/base/port/labgrs.html">LabGRS</a> | Noviembre 2023]

---
class: center,middle
background-image: url(data:image/png;base64,#labgrs_logo.png)
background-size: 35%

---
## Contenidos

.pull-left[
- Combinaciones de bandas.

- Calibración radiométrica. (TOA radianza - TOA reflectancia)

- Corrección atmosférica.

- Filtros de calidad de la información satelital.

- Compensación topográfica por iluminación.
 
 ]
 
.pull-right[
<img src="data:image/png;base64,#https://raw.githubusercontent.com/allisonhorst/stats-illustrations/main/rstats-artwork/r_rollercoaster.png" width="650px"/>
© Allison Horst
]

---
### Recordemos.

<center>
<img src="data:image/png;base64,#Signatura_vegetacion.PNG" alt="drawing" width="800" height="480"/>
<center>
---

---
### 1) Combinaciones de bandas.

A partir de las signaturas espectrales es posible entender el comportamiento de los objetos frente a diferentes longitudes de onda.

Esto nos permite crear combinaciones de bandas en los canales **RGB** para resaltar elementos específicos en superficie.

```r
#librerias necesarias ---- 
library(satellite)
library(tidyverse)
library(terra)
library(raster)
library(landsat)
list_landsat <- list.files(path = "L08_Concepcion_C02_LVL1/", 
 pattern = glob2rx("*T1_B*.TIF"),
 full.names = T) # lista de bandas Landsat 8 dentro del directorio.
list_landsat_plot <- list_landsat[-c(2, 3, 10, 11)] # trabajar solo con las bandas 1-7.
l8_st <- list_landsat_plot %>% rast() # crear un raster solo con bandas a 30 m.

plotRGB(l8_st, r=4, g=3, b=2, stretch = "lin", axes = TRUE)
plotRGB(l8_st, r=4, g=5, b=2, stretch = "lin", axes = TRUE)
plotRGB(l8_st, r=5, g=6, b=7, stretch = "lin", axes = TRUE)
```

---
### 1) Combinaciones de bandas.

<center>
<img src="data:image/png;base64,#comb_bandas.PNG" alt="drawing" width="1000" height="500"/>
<center>
---
### Abrir imagen satelital con la librería satellite.

> "El objetivo principal de la librería es proporcionar clases y métodos estándar para agilizar el flujo de trabajo de sensores remotos en R. Proporciona su propia clase de objeto (satélite) junto con métodos estándar para transformaciones básicas de imágenes, como correcciones atmosféricas y topográficas, entre otras." (Nauss & Detsch, 2021)

para mas informacion visitar la descripción del paquete [**satellite**](https://CRAN.R-project.org/package=satellite)

```r
list_landsat <- list.files(path = "L08_Concepcion_C02_LVL1/", 
 pattern = glob2rx("*T1_B*.TIF"),
 full.names = T) # lista de bandas Landsat 8 dentro del directorio 
l8 <- satellite(list_landsat) # objeto satellite / capas, metadato, bitacora 
l8 <- l8[[1:7]] # trabajar solo con las bandas 1-7 
class(l8) # vemos el tipo de objeto. 
```

```
## [1] "Satellite"
## attr(,"package")
## [1] "satellite"
```

---
### Objeto tipo “satellite”.

Un objeto de tipo satellite se crea a partir de una lista de archivos asociados al directorio de mi información satelital. Estos objetos cuentan con 3 elementos internos:

1. **Layer**: objeto de tipo ráster layer el cual almacena las bandas de mi satélite.  
 2. **Meta**: objeto de tipo data.frame que almacena la información contenida en nuestro metadato. (Archivo MTL.txt)

3. **Log**: objeto de tipo lista que almacena los cambios de nuestro objeto satellite. (correcciones radiométricas, atmosféricas etc.)

---
### 2) Calibración radiométrica.

---
### 2) Calibración radiométrica.

En un gran número de sensores (TM, ETM+, OLI, etc.) se establece la conversión a **radianza espectral** en base a una **relación lineal** existente entre el nivel digital bruto y los valores de radianza (unidades físicas). Sin embargo, **siempre es importante** revisar los manuales de usuario de cada producto con el que estemos trabajando.

De tal forma que la ecuación que describe la radianza espectral, según Chander et al. (2009) es:

---
### 2) Calibración radiométrica

```r
## ND a TOA radianza calculo manual -----
metadato <- l8@meta # Crear tabla con los metadatos 
meta_gain_offset <- data.frame(gain = metadato$RADM, offset = metadato$RADA, banda = metadato$TYPE)
print(meta_gain_offset)
```

```
##         gain    offset banda
## 1 0.01293900 -64.69457   VIS
## 2 0.01325000 -66.24803   VIS
## 3 0.01220900 -61.04700   VIS
## 4 0.01029600 -51.47827   VIS
## 5 0.00630040 -31.50213   NIR
## 6 0.00156690  -7.83429  SWIR
## 7 0.00052812  -2.64058  SWIR
```

1. **gain**: factor multiplicativo.  
2. **offset**: factor aditivo.

---
### Niveles digitales a radianza (TOA).

```r
st <- l8@layers %>% stack # la compatibilidad de los objetos satellite es con raster.
rad1 <- st[[1:7]]*meta_gain_offset$gain[1:7]+meta_gain_offset$offset[1:7]
```

Calibración radiométrica. Funcion convSC2Rad (basado en los manuales del USGS)

```r
rad2 <- convSC2Rad(l8) 
```

**Descripción:** Convierte los niveles digitales brutos de mi imagen en **radianza** usando una conversión linear simple sin ningún tipo de corrección atmosferica aplicada.

---
### Niveles digitales a radianza (TOA).

---

<center>
<img src="data:image/png;base64,#ND_TOA_REF.PNG" alt="drawing" width="1050" height="600"/>
<center>
---
### Niveles digitales a reflectancia aparente (TOA reflectancia)

- Transformar la radianza a reflectancia al tope de la atmosfera es un paso intermedio a la corrección atmosférica final.

- Si bien la mayoría de los métodos de corrección transforman automáticamente los niveles digitales a reflectancia superficial, es importante conocer el **porcentaje** de energía que es reflejada por los objetos en la superficie y que finalmente es captada por el sensor.

- Al ser la reflectancia al tope de la atmosfera, se le conoce como **reflectancia aparente**.
---
<center>
<img src="data:image/png;base64,#formulas_TOA_REF.PNG" alt="drawing" width="1100" height="630"/>
<center>

---
### Niveles digitales a reflectancia aparente (TOA reflectancia).

```r
meta_mult_add<- data.frame(add = metadato$REFA, mult = metadato$REFM, banda = metadato$TYPE)
print(meta_mult_add)
```

```
##    add  mult banda
## 1 -0.1 2e-05   VIS
## 2 -0.1 2e-05   VIS
## 3 -0.1 2e-05   VIS
## 4 -0.1 2e-05   VIS
## 5 -0.1 2e-05   NIR
## 6 -0.1 2e-05  SWIR
## 7 -0.1 2e-05  SWIR
```

```r
sun_elev <- metadato$SELV[1] # extraer angulo elevación solar del metadato.

st <- l8@layers %>% stack # transformar mi objeto satellite a stack. 
TOA_ref_zenith <- (st[[1:7]]*meta_mult_add$mult[1:7]+meta_mult_add$add[1:7])/
 (sin(sun_elev*3.1416/180)) # ND a reflectancia aparente. 
```
---
### Niveles digitales a reflectancia aparente (TOA reflectancia).

```r
TOA_ref_dn <- convSC2Ref(l8, szen_correction = "TRUE")
```

- De esta manera la función automáticamente detecta los parámetros internos en los metadatos.

- En caso de utilizar **szen_correction = F**, no se realizará la corrección por ángulo solar.

---

```r
par(mfrow = c(1,3))
plot(l8@layers[[2]], main = "niveles digitales ND") # BANDA AZUL
plot(rad2@layers[[13]], main = "Radianza TOA") # BANDA AZUL
plot(TOA_ref_dn@layers[[13]], main = "Reflectancia TOA top of atmosphere") # BANDA AZUL
```

---

### 3) Corrección atmosferica.

Clasificación de correcciones atmosféricas:

- **Mediciones in situ**

- **Mixtas** Ejemplo: Landsat Ecosystem Disturbance Adaptative Processing System (LEPADS) y Land Surface Reflectance Code (LaSRC)

- **A partir de modelos físicos de transferencia radiativa** Ejemplos: Second Simulation of a Satellite Signal in the Solar Spectrum (6S vector version), Multi-angle implementation of atmospheri correction (MAIAC).

- **Con datos de la propia imagen** Ejemplos: Simple Dark object Substraction (sDOS), Sen2cor, ACOLITE.
---
<center>
<img src="data:image/png;base64,#reF_app_SR.PNG" alt="drawing" width="1000" height="600"/>
<center>

---
<center>
<img src="data:image/png;base64,#softwares.PNG" alt="drawing" width="900" height="500"/>
<center>

---
### 3) Corrección atmosferica. 
 Métodos darkest object substraction:

En la literatura podemos encontrar diferentes métodos DOS los cuales se basan en principios similares.

El paradigma inicial señala que algunos pixeles son cuerpos oscuros cuya reflectancia seria cercana a 0. Esto indica que la radianza captada por el sensor de estos pixeles oscuros, estaría influida por el efecto de la **bruma**.

Para estos métodos se estima que los piexeles oscuros reflejan al menos entre el 1-2% de la energía. (0.01)

El método **darkest pixel subtraction** busca reducir o eliminar este efecto de la bruma para obtener la reflectancia superficial.
---

---
<center>
<img src="data:image/png;base64,#METODOS.PNG" alt="drawing" width="730" height="650"/>
<center>
---
### Metodo DOS2 Chavez 1996.

Este método de corrección por DOS presento una mejora respecto al **simple darkest object subtraction** al tener en consideración la corrección de la transmitancia atmosférica desde la superficie terrestre hacia el sol.

También conocido como método CosTz, esta mejora utiliza el ángulo zenital del sol el cual está asociado a la transmitancia atmosférica para las fechas especificas en las cuales se tomó la imagen.

Sin embargo este método sigue considerando la Transmitancia de la energía del sol a la tierra (Tv) como constante (no pierde energía en el viaje de entrada se mantiene con valor 1 como el caso del sDOS).

---

---
### Metodo DOS2 Chavez 1996.

```r
# reflectancia superficial metodo DOS1:
# Lp = efeccto bruma / path radiance / Lhaze
# Tv = transmitancia de la atmosfera direccion vision
# Tz = transmitancia atm en direccion de la iluminacion
# Edown = irradiancia difusa descendiente

# Lp = Ldark - Ldo1% 
## Ldark = G * DNmin + B
```

```r
df_st <- l8_st %>% as.data.frame() %>% na.omit() # transformar mi stack a data frame para obtener el valor de pixel mas oscuro
DNmin <- c(8731, 7720, 6293, 5500, 4040, 4067, 4795) # pixeles oscuros para cada banda
# calcDODN() función del paquete satellite para obtener los ND mas oscuros por banda
G <- metadato$RADA # factor aditiv
B <- metadato$RADM # factor multi
Ldark = B[1:7]*DNmin[1:7]+G[1:7] # valores minimos por banda
```

---
Calcular la irradianza del sol a partir de las radianzas y reflectancias aparentes maximas.

```r
## Esun = (pi * d2) * RadMax / RefMax         / esto para cada banda. 
d = metadato$ESD[1] # distancia al sol [AU]
sun_elev = metadato$SELV[1] # angulo elev sol
RadMax = metadato$RADMAX # Radianza maxima
RefMax = metadato$REFMAX # Reflectancia maxima
zenit = metadata$SZEN[[1]] # angulo zenit del sol

Esun = (3.1416*d^2)*RadMax[1:7]/RefMax[1:7]
```

```
## [1] 1972.25821 2019.61650 1861.05925 1569.35010  960.36396  238.83378   80.49978
```

Estimar la radianza del sensor.

```r
### ND - radianza (TOA)
meta_gain_offset <- data.frame(gain = metadato$RADM, offset = metadato$RADA, 
 banda = metadato$TYPE)
print(meta_gain_offset)

rad1 <- l8_st[[1:7]]*meta_gain_offset$gain[1:7]+meta_gain_offset$offset[1:7]
```
 
---
Estimar Lp / Efecto de la bruma sobre cada banda.

```r
Edown = 0
Tv = 1
Tz = cos(zenit*3.1416/180)

Lp <- Lmin[1:7]-0.01*(Esun[1:7]*cos(zenit*3.14/180)*Tz+Edown)*Tv/(3.1416*d^2)
```

```
## [1] 43.3994666 31.0485051 11.1828020  1.2695405 -8.4229964 -2.0522199 -0.3072788
```

Correccion atmosferica metodo DOS1

```r
SR <- 3.1416*(rad1[[1:7]]-Lp[1:7])*d^2/(Tv*Esun[1:7]*cos(zenit*3.1416/180)*Tz)+Edown
```

---
<center>
<img src="data:image/png;base64,#ND_SR.PNG" alt="drawing" width="1000" height="600"/>
<center>

---
### Metodo DOS2 Chavez 1996.

A partir de la función **calcAtmosCorr** podemos realizar el método de DOS2 el cual sigue la misma filosofía de la ecuación de reflectancia inicial con un cambio en los parámetros.

```r
list_landsat <- list.files(path = "L08_Concepcion_C02_LVL1/", 
 pattern = glob2rx("*T1_B*.TIF"),
 full.names = T) 
l8 <- satellite(list_landsat) # objeto satellite / capas, metadato, bitacora 
l8 <- l8[[1:7]]
DOS2_SR <- calcAtmosCorr(l8, model = "DOS2", esun_method = "RadRef")
```

Para el cálculo de la irradiancia solar (Esun) se utiliza el mismo método aplicado en el DOS1 considerando la radianza máxima y reflectancia máxima de cada banda.

![A GitHub Octocat, style "mother hubber"](data:image/png;base64,#Esun.png)
---
### Bibliografía.

Thomas Nauss, F.D. (2021) satellite: Handling and Manipulating Remote Sensing Data, Satellite - classes and methods for satellite data.

Chavez, Jr, Pat. (1996). Image-Based Atmospheric Corrections - Revisited and Improved. Photogrammetric Engineering and Remote Sensing. 62. 1025-1036.

Chavez & P.S, 2015. Image-Based Atmospheric Corrections - Revisited and Improved. , no. January 1996.

Moran, M.S., Jackson, R.D., Slater, P.N. y Teiuet, P.M., 1992. Evaluation of Simplified Procedures for Retrieval of Land Surface Reflectance Factors from Satellite Sensor Output. , vol. 184, pp. 169-184.

Song, C., Woodcock, C.E., Seto, K.C., Lenney, M.P. y Macomber, S.A., 2000. Classification and Change Detection Using Landsat TM Data : When and How to Correct Atmospheric Effects ? , vol. 4257, no. 00.

Wang, Y. , Wang, X. y He, H., 2019. An Improved Dark Object Subtraction Method for Atmospheric Correction of Remote Sensing Images [en línea]. S.l.: Springer Singapore. ISBN 9789811399176. Disponible en: http://dx.doi.org/10.1007/978-981-13-9917-6_41.

---
class: inverse middle
![logo](data:image/png;base64,#LABGRS-PUCV.png)